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Abstract 


A novel Computational Fluid Dynamics (CFD) coupling framework that uses a conventional 
Reynolds-Averaged Navier-Stokes (RANS) solver to resolve the flow field near the body and 
Particle-based Vorticity Transport Method (PVTM) to predict the evolution of the far field 
wake is developed, refined, and evaluated for fixed and rotary wing cases. For rotary wing 
case the RANS/PVTM modules are loosely coupled to a Computational Structural Dynamics 
(CSD) module that provides blade motion and vehicle trim information. The results from the 
coupled framework are compared with several experimental data sets (a fixed-wing wind 
tunnel test and a rotary-wing hover test). The PVTM module is refined by the additions of 
vortex diffusion, stretching, and reorientation models as well as an efficient memory model. 
Validation with the fixed-wing wind tunnel test data shows that the coupled RANS/PVTM 
method provides good prediction on wing performance (pressure distribution and sectional 
loads) and tip vortex parameters (core size, location and swirl velocity). For the rotary wing 
under hover condition, the results from the coupled RANS/PVTM/CSD framework correlate 
well with the hover test data, however the simulation over-predicts the rotor torque by about 
50%. Overall, the tip vortex parameter validations are good. The tip vortex swirl velocity is 
slightly over-predicted, and the difference in tip vortex trajectory is within one chord length 
over 150° wake age. Significant improvement on the correlations of vortex trajectory and 
swirl velocity is observed when the vortex stretching model is used in the PVTM module. 
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Chapter I 


Introduction 


1.1 Background and motivation 

Accurately predicting the long-term dynamics of the wake produced by wings and rotors under 
high Reynolds number flow condition is still one of the most challenging tasks for CFD 
simulations. The wake structure behind a wing under such condition, in the simplest form, 
consists of vortex sheet from the boundary layer and a strong tip vortex from 3D finite wing 
effects. The tip vortex is quickly formed and is normally well organized before leaving the wing 
trailing edge [1.1]. Typical velocity profiles of the wing tip vortices can be found in Ref. [1.2]. 
The tip vortex continues strengthening after leaving the trailing edge by rolling in the vortex 
sheet. Then the viscous diffusion effect causes the tip vortex to grow in size and lose strength 
very slowly as it convects downstream. 

The wake structure from a rotor is much more complex since all the vortex sheets, tip, and root 
vortices released from all of the blades are intertwined due to the rotation of the blades to form a 
helical wake structure [1.3]-[1.4]. Unlike the wake from wings, the helical structure of rotor 
wake produces very strong blade-wake and wake-wake aerodynamic interactions. The blade- 
wake aerodynamic interaction often causes the unpleasant Blade Vortex Interaction (BVI) noise 

[1-5]. 

Many researchers have simulated the wing and rotor wake flow successfully using a variety of 
classical and physics-based approaches. Simulations of the wake flow behind wings were 
attempted using vortex filament method [1.6], vortex panel method [1.7], a classical wing theory 
[1.8], Computational Fluid Dynamics (CFD) [1.9], and vortex particle method [1.10]. The CFD 
calculations in Ref. [1.9] used two CFD solvers (Reynolds-Averaged Navier-Stokes (RANS) and 
hybrid Euler/Large Eddy Simulation (LES) solvers) to compute aircraft wake very far 
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downstream, with an equivalent distance of more than two nautical miles (~500c) for a typical 
commercial aircraft. With the vortex particle method [1.10], a high resolution wake from an 
aircraft was simulated using billions of vortex particles. 

Four major approaches were used to simulate rotor wake flow - (i) vortex lattice, (ii) vortex blob 
or particle, (Hi) full Computational Fluid Dynamics (CFD), and (iv) Vorticity Transport Method 
(VTM). The first and most widely used method for rotor wake prediction is using vortex lattices 
to represent the shed and trailed vorticity generated by the rotor [1.3]-[1.4]. Both straight and 
curved vortex lattices have been considered. Generally, these vortex elements are free to move 
with local velocity induced by other vortex elements in the flow field, however less sophisticated 
rigid wake models use empirical data to define the location of the lattices and do not allow the 
vortex elements to move. In either case the induced velocity is calculated using the Biot-Savart 
law, which is a computationally expensive operation. In order to minimize the amount of 
computation, the trailing vortices are often approximated as a single or few tip vortices after 60° 
to 90° wake age and the shed vortices are kept up to 30° to 45° behind the rotor. The strength of 
the vortex lattices is calculated from the bound circulation on the blade, which is derived from 
the blade sectional lift. Other parameters for this vortex lattice method such as core models and 
releasing points are very difficult to determine. 

Vortex particle method was also used to capture shed and trailed vorticity in the rotor wake 
[1.1 1]-[1.12]. Sometime the term vortex “blob” is introduced as collections of vortex particles. 
The vortex particles or blobs can also move with local velocity which is calculated again using 
the Biot-Savart law for every particle in the computational domain. Since a large number of 
particles (N) is needed to define the entire flow field the computational cost can be as high as N 
operations although some techniques can be used to reduce this cost, e.g. Particle-Mesh (NlogN) 
and Fast Multipole (NlogN). The strength of vortex particles again can be determined from the 
bound circulation on the blade. Core models of this vortex particle method are well established, 
but the location of releasing points is still challenging to obtain. 

Full CFD has also been considered for simulating the rotor wake [ 1 . 13]-[ 1.14]. This approach 
solves the Navier-Stokes equations for primitive variables (i.e. velocity, density, and pressure) in 
the flow field at all grid points in the entire domain. Normally, two sub domains are 
implemented in rotorcraft CFD modeling. One domain rotates with the blades and another non- 
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rotating domain is used in the background to capture the rotor wake. Interpolation between the 
two domains has to be performed at every time step. The rotating domain usually extends three 
to five chord lengths from the rotor surface, while the background grid extends three to five rotor 
radii from the surface. This method has the fidelity to capture boundary layer and 3-D effects 
near the root and tip of the blade very well. However, a very large number of grid points is 
necessary to simulate the vortex wake satisfactorily, since the cell dimensions inside the 
boundary layer can be as small as 1/10000 of the chord length and the cell dimensions in the 
vicinity of the tip vortices can be approximately 1/50 of the chord length. Grid adaptation is a 
necessity to properly distribute more grid points to the regions with high vorticity, and these 
regions usually move constantly throughout the field. 

A rotor wake can also be modeled using VTM [1.15], The vorticity transport equations, which 
are derived from conservation of momentum, are solved to determine the evolution of vortex 
parameters on a Cartesian grid. The maximum resolution of the grid typically used for this 
approach is about a quarter of the chord length. Similar to other vortex approaches, the induced 
velocity is calculated using the Biot-Savart law. A multi-level grid is introduced to reduce the 
computational cost. The initial strength of the vortex is obtained from 2D lifting line theory. 
This VTM approach can satisfactorily model the evolution of the rotor wake. 

Overall traditional RANS solvers can predict the generation of vorticity very well, but often have 
difficulties in simulating the evolution of the vorticity field over a long period of time because of 
limitations due to numerical dissipation and regeneration/adaptation of the grid. On the other 
hand, the evolution of the vorticity field can be modeled satisfactorily using VTM or modem 
vortex particle method which are both based on the vorticity transport equations. In addition, the 
vortex particle method does not require grid adaptation, since there is a natural tendency of 
particles to cluster around regions with high vorticity. Therefore, a hybrid approach using a 
RANS solver near the surface to calculate the generation of vorticity and a Particle-based VTM 
(PVTM) to simulate the evolution of the far field wake flow offers a unique capability to capture 
the entire wing or rotor wake flow with high accuracy. 

Recently, a new hybrid approach using fully coupled Reynolds-Averaged Navier-Stokes (RANS) 
and Particle-based Vorticity Transport Method (PVTM) solvers was introduced to simulate 
rotorcraft wake flow [1.16]. The approach divides the flow field into several regions and uses 
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appropriate flow solvers according to the dominant physical features of the flow in each region. 
The first region covers the flow field near aerodynamic surfaces (extending about one chord 
from the surfaces), where the flow features are dominated by the effects of boundary layer 
viscosity and the geometry of the airfoil and blade planform. The near body flow field is 
resolved using a 3D compressible RANS solver [1.17]. Outside of the RANS regions, the flow 
field is primarily dominated by the vortices being shed from the aerodynamic surfaces. This 
vortex-dominated flow region is simulated using a Particle-based Vorticity Transport Method 
(PVTM). The influence of this far field flow region is transferred to the RANS regions using the 
field velocity approach [1.17]. By modeling the far field with PVTM, the shed vorticity can 
remain well organized and in particular the tip vortices can maintain their compact and stable 
cores for a longer period of time than vortices simulated using conventional RANS/CFD. This 
approach is particularly suitable for rotorcraft applications where maintaining the correct 
velocity gradient in the vortex core is critical for acoustic and vibratory loads calculations, which 
depend greatly on the location, size, and strength of the vortex. A previous study showed that 
the methodology is viable, but did not illustrate quantitatively how well the approach predicted 
the near body flow field and associated vortex-dominated flow field [1.16]. 

1.2 Focus of the present research 

Further development and validations of the fully coupled RANS/PVTM/CSD methodology 
are achieved in this study (Fig. 1.1). The development efforts include - (i) refining coupling 
methodology between RANS/PVTM, (ii) coupling the RANS/PVTM with Computational 
Structural Dynamics (CSD) code, (iii) implementing a diffusion model for PVTM, (iv) adding a 
vortex stretching model for PVTM, (v) implementing a dynamic memory model for PVTM to 
reduce memory usage. The coupled RANS/PVTM methodology is validated against several 
experimental data sets including - (a) a low aspect ratio fixed-wing wind tunnel test, and (b) a 
hover test for a lightly loaded rotor. 
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Chapter II 


Coupled PVTM/RANS/CSD Methodology 


2.1 Coupled RANS/PVTM methodology 

In the present study, the entire flow field is divided into near and far fields, which are resolved 
using RANS and PVTM solvers, respectively. The RANS solver predicts the generation of the 
vorticity field near aerodynamic surfaces and releases the vorticity into the PVTM domain 
through a convection process at the boundary of the RANS domains. The PVTM solver 
determines the evolution of the entire far field vorticity after the vortex particles are released 
from the RANS domains. At every time step, the RANS domains provide vorticity information 
to PVTM, and PVTM provides current induced velocity information from all vortex particles in 
the far field to the near field RANS grids. A summary of the RANS and PVTM methodologies 
and the coupling process is provided in the following sections. 

2.1.1 Near field analysis (RANS domain) 

A RANS solver is used to calculate the generation of the vorticity field near the aerodynamic 
surfaces (wing or rotor blades), and in the present study the solver is based on a 3D RANS 
implementation described in Ref. [1.17]. The analysis solves the RANS Equations for primitive 
variables (i.e. p , V, and P ) inside a 3D structured grid using a finite volume approach. This 
RANS domain defined by the grid usually extends to approximately one chord length away from 
the surface. The short distance to the boundary minimizes numerical dissipation of the generated 
vorticity field before reaching the boundary and releasing into the PVTM domain. The boundary 
conditions for this RANS domain are - (i) the induced velocity from the far field (PVTM and 
other disconnected RANS domains), (ii) the free stream density, and (iii) the non-reflective 
boundary condition (extrapolation of pressure from inside the RANS grid). In addition, the 
induced velocity from the far field (PVTM and other disconnected RANS domains) is included 
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as field velocity [1.17] to all of the grid points in the RANS domain to ensure accuracy and to 
speed up convergence of the coupling process. 

For the fixed wing computation, only one RANS domain is used. The RANS grid extends to 
about one chord in all directions except behind the trailing edge, where it extends only to about a 
half chord length to minimize the dissipation of vorticity to be transferred into the PVTM 
domain. An example of the near body RANS grid for a fixed wing calculation is shown in Fig. 
2.1a. The RANS grid has a dimension of 249x131x79 grid points. 


The near field calculation with a RANS solver in the rotor simulation involves as many RANS 
domains as the number of blades in the rotor. All domains are disconnected from the other 
RANS domains. Figure 2.1b shows the RANS grid system required for a four bladed rotor 
simulation. The coverage of the grid points is similar to the RANS grid for the fixed wing case: 
0.5c behind the trailing edge, and lc in other directions. For each blade, the RANS grid has a 
dimension of 239x151x69 grid points. The influence of one blade to the other blades is 
calculated by - (i) converting the velocity field inside the RANS domain into a vorticity field, (ii) 
converting the vorticity field into a particle vortex field, and (iii) calculating the additional 
induced velocity from this near blade particle vortex field to other blades. 

2.1.2 Far field analysis (PVTM domain) 


Outside of the near surface RANS domains, the flow field is represented by collections of three- 
dimensional vortex particles similar to those presented in Ref. [2.1]. Each vortex particle has 
two vector quantities, location and strength, associated with it. The strength of the vortex is a 
volume integration of the vorticity field around the particle, and the evolution of the strength is 
governed by the vorticity transport equations: 


J ^ dV t = J [a> ■ Vu]dV, + J [vV 2 ©] cIV, + S[ 

a=\co dV t 


— = a- Vm + a H + a,, 
dt d ' 

dx 

— = u 

dt 


( 2 . 1 ) 

( 2 . 2 ) 

(2.3) 
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where u is the local flow velocity, 0) is vorticity, v is kinematic viscosity, and S , is the vorticity 
released from the RANS domain into the differential volume dVj. Using the vortex particle 
strength definition in Eq. (2.2), Eq. (2.1) becomes the vortex particle governing equation (Eq. 
2.3), where oris the strength of a vortex particle in dVi, a‘Vu is the vortex stretching term which 
cannot incur reduction or increase in the vorticity strength, |«| [2.2], a d is the rate of change of 
tge vortex particle strength due to diffusion, a r is the rate of release of vortex particle strength 
from RANS domains, and x is the location of the vortex particle. The differential volume, dV h 
represents the PVTM cell resolution, and a 3D Cartesian volume is chosen for convenience. 
Representations of the multi-level PVTM computational domains are presented in Fig. 2.2 for 
the fixed wing and the rotor cases. To reduce computational cost, only one vortex particle is 
allowed in a PVTM cell. When two or more particles are inside the same PVTM computational 
cell, they are merged into one particle as graphically shown in Fig. 2.3b, with the new strength 
being the sum of the strength of all particles, and the location being the strength weighted 
centroid of all the particles. The location of each vortex particle is governed by the local velocity 
induced from - (i) all vortex particles in the far field, (ii) vortex particles representing vorticity in 
the RANS domains. The detailed derivation of vortex particle evolution equations are given in 
Chapter III, while the derivation of the released vortex from RANS domains is given below. 

2.1.3 Coupling RANS/PVTM domains 

In order to couple the RANS and PVTM domains, proper information is passed between them at 
every time step. The information passing from RANS domains to the PVTM domain is the 
vorticity field for calculating - (i) the vortex particles released from RANS domains to PVTM 
domain, (ii) the induced velocity from RANS domains. The information passing from PVTM 
domain to RANS domains is the vortex particle field which represents the entire far field 
vorticity that is used to calculate the far field induced velocity. Additional information transfer 
(vortex particles representation of the RANS vorticity field) between the RANS domains is 
needed for calculating the induced velocity from one RANS domain to the other RANS domains. 

2. 1.3.1 Vortex particles released from RANS domains to PVTM domain 

The vorticity released from RANS domains into the PVTM domain, a,-, is derived from the 
convection of the vorticity field from the boundary of the RANS domains. For each RANS cell 
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that is on the boundary, the vorticity convected from the RANS domain into the PVTM domain 
is defined as follows: 


a r = • (n ■ n) dt (2.4) 

X r= X mid+\ U dt 

where a, is the strength of the released particle, X r is the location of the released particle, A is 
area of the outer boundary, ft is the outward normal vector and X mi( j is the mid-point of the area 
A. In addition to these convectively released vortex particles, the moving RANS grid also 
releases vortex particles due to the movement of the grid boundary as follows: 

a, = j <Mu dt = '/{at + af~')A \Ku - K{\ (2-5) 

x r =y,(xL, + KA 

After the release, the vortex particles in the same PVTM computational cell are merged to reduce 
the computational cost. Graphical representation of the vorticity releasing process from RANS 
domain to PVTM domain is presented in Fig. 2.4. 

2. 1.3.2 Induced velocity from RANS domains to PVTM domain 

The aerodynamic influence from RANS domains to the PVTM domain is taken into account by 
including the induced velocity from vorticity in the RANS domains to calculate the particle 
velocity in Eq. (2.3). The continuous vorticity field in the RANS domains is processed into a 
particle vortex field representing the RANS domains. These vortex particles representing the 
RANS domains are included in calculation of the local flow velocity for vortex evolution 
equations. Detailed derivation of the induced velocity calculation is given in Chapter III. 

2. 1.3.3 Induced velocity from PVTM domain to RANS domains 

The far field induced velocity from the PVTM domain is included to the near field RANS 
calculation as field velocity [1.17]. The induced velocity is added to all of the grid points in the 
RANS domain to ensure accuracy and to speed up convergence of the coupling process. To 
reduce the computational effort, the induced velocity from the PVTM domain is calculated on an 
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interpolation grid and the induced velocity on all RANS grid points are calculated from 3D linear 
interpolation from the induced velocities on the interpolation grid (see Fig. 2.5) . 

2. 1.3.4 Induced velocity from a RANS domain to other RANS domains 

When more than one RANS domains (not connected to each other) are used, the influence 
between these domains are accounted for by calculating the induced velocity from one RANS 
domain to the other RANS domains. Again, the continuous vorticity field in a RANS domain is 
converted to a vortex particle field that represents the vorticity in the RANS domain. The 
induced velocity from this vortex particle field is included in the calculation of the induced 
velocity for other RANS domains. 

2.2 Coupling RANS/PVTM with CSD 

For the rotor case, a loose coupled trim methodology is used to couple a CSD code (CAMRAD 
II [2.3]) with the RANS/PVTM analysis. The comprehensive CSD code, CAMRAD II, handles 
the vehicle trim calculation and performs a structural analysis for the blade motion in response to 
the RANS/PVTM aerodynamic loading. The loose coupled trim methodology presented in Ref. 
[2.4] is adopted for coupling CAMRAD II and RANS/PVTM and is summarized in Fig. 2.6. 

The calculations for the vehicle trim and the structural analyses are performed by CAMRAD II 
with the aerodynamic loading from RANS/PVTM analysis in a once-per-revolution basis (a 
more frequent information transfer is also possible). The pressure and shear forces on the blades 
(from the near blade RANS domains) are integrated to provide the blade sectional load ( M 2 Cl , 
M 2 C d , and M 2 C m ) at various collocation points along the blade quarter chord line for the CSD 
analysis. The CSD code uses this blade sectional loading to retrim the rotor and provides the 
RANS/PVTM analysis a new blade motion (6 DOFs) at the collocation points. The new blade 
motion is used in the RANS/PVTM analysis to calculate the new blade aerodynamic loadings. 
This coupled trim process continues until the following variables are converged: (i) trim 
parameters, (ii) blade motion, and (iii) blade loading. 
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Chapter III 


PVTM Methodology 


3.1 Vorticity Transport Equations and solution methodology 

Vorticity transport equations (derived from taking the curl of the momentum equations from 
Navier-Stokes Eqs.) can be expressed in several frames of reference - (i) Eulerian frame [1.15] or 
(ii) Lagragian frame [2.1]. In the current study, the Lagragian frame implementation is adopted 
for convenience since it offers gridless capability (no flow calculations at the grid points). The 
approach discretized the vorticity field into a vortex particle field, and each vortex particle 
carries a local vorticity field around it. The movement of these vortex particles represents the 
transportation of vorticity without any loss in vorticity strength. Again, the vorticity transport 
Eq. (2.3) is presented here in a slightly different form: 

a k+ ' = a k + jV • Vw) dt + a d + a r (3. la) 

x M =x k +\udt (3.1b) 

where superscript k represents a time step index, a= foodV, u is the local velocity, is the 
change in vortex particle strength due to the diffusion effect, a,- is the strength of the vortex 
particle released from the RANS domains (Eqs. 2.4-2. 5), and x is the location of the vortex 
particle. Equation 3.1a governs the change in the strength of the vortex particle in a particular 
PVTM cell, while the location of the particle is calculated using Eq. 3.1b. The detailed 
definition of each term in Eq. 3. 1 can be found in Section 2.1.2. 

Equations 3.1 are solved using explicit integration schemes based on Runge-Kutta (RK) 
integration method with 4 th or 6 th order accuracy in time [3.1]. Equation 3.1b is solved directly 
with the RK scheme with adaptive time step (see Appendix A). When solving Eq. 3.1a, the 
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terms a & and a,- are calculated in separate steps and included at the end of the time step, thus 
without these terms, Eq. 3.1a reduces to: 

a k+1 = (I + A e )a k (3.2) 

A £ = | Vw dr 

where A e is the strain tensor of the local velocity field. It should be noted that the summation of 
the eigenvalues of A e is zero (continuity condition). There is also one eigenvalue of A e that is 
zero, Ao, representing a solution that conserves the total vorticity (\ct +1 \ = \ct\). Other 
eigenvalues are associated with non-physical increase or decrease in vorticity strength. Thus, the 
only solution that conserves vorticity must have the strength vector oriented in the direction of 
the eigenvector associated with A 0 . Numerically, the solution of Eq. 3.2 represents a 
reorientation or 3D-rotation of the vortex strength vector, ot, into the direction of the eigenvector 
associated with Ao. The derivation of this eigenvector is given in Appendix B for a system of 
two vortex particles. Due to the fact that the strain tensor is a summation of the stain field 
induced from all vortex particles in the entire field, a superposition method can be used to 
determine the reorientation of the strength vector of a vortex particle in response to other vortex 
particles in the field. 

3.2 Induced velocity calculation and computational scheme 

Calculation of induced velocity involves summation of the induced velocity from all vortex 
particles in the flow field (far field PVTM domain and near field RANS domains). The induced 
velocity at an arbitrary location X from N vortex particles can be calculated using the modified 
Biot-Savart law [2.1]: 


u(X) = ± 


i = 1 


1 (X-X^XOC, 

4 *{x-x,\ 2 +S 2 f 


(33) 


where 8 is the desingularized parameter which is introduced to avoid singularity of the induce 
velocity when X = x,. The value of 8 is set to be 0. lxaf. The numerical effect of this 
desingularized parameter is similar to using a vortex core radius of 8. 
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The calculation of the induced velocity is very computationally expensive. In order to reduce the 
number of induced velocity calculations within the RK integration substep, the induced velocity 
is calculated on tetrahedral basis points (Fig. 3.1) and a trilinear interpolation is used to 
determine the induced velocity within the RK substep. The distance from the particle to the all 
vertices of the tetrahedral is set to 0.5% of the PVTM cell size. The interpolation is based on the 
first order Taylor expansion of the quantity u about the center of the tetrahedral. The derivatives 
of the quantity u are calculated using Eq. (3.4) 


du 

dx 

du 

dy 

du 

dz 



(3.4) 


where the subscript i represents the i th vertex of the tetrahedral, and the subscript o represents the 
values at the center of the tetrahedral ( u 0 is the average value of uj, ui, uj, 114 ). The interpolated 
value, u'(x'y'z% is obtained following Eq. (3.5) 


u =u a + 





x'-x a 

dll 

du 

du 

0 

/ 

dx 


dz 

y -y Q 




L z ~ z o] 


(3.5) 


3.3 Local strain tensor calculation 


The calculation of the local strain tensor is straightforward since the strain rate tensor is already 
calculated from Eq. (3.4) for the velocity interpolation. A first order time integration scheme is 
used to determine the strain tensor in Eq. (3.2b). 

3.4 Reorientation of vortex strength vector 

The strength vector of the vortex particles is reoriented or rotated toward the direction of the 
eigenvector associated with Ao for a vorticity solution that conserves the total vorticity. It is 
shown in Appendix B that the reorientation of the vortex strength vector can be achieved by 
using a vortex stretching term in the following form: 
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(3.6) 


a 


4 K 


(a r x«r i )x«r / ’ 


where ij dictates how fast the reorientation of the strength vector occurs. An adjustment in the 
velocity calculation is required to mitigate the numerical error accumulation due to the explicit 
time integration scheme (see Appendix B). 

3.5 Diffusion model for PVTM 


3.5.1 Diffusion model using vortex redistribution 

For vortex method, the effect of diffusion process can be simulated using a vortex redistribution 
method. In this study, the diffusion effect is simulated by splitting a vortex particle into five 
smaller particles and redistributing the original vortex strength to the new particles. The 
locations of these five particles are defined by the vertices of a randomly oriented tetrahedral and 
its center (see Fig. 3.2). The distance from the original particle (center of the tetrahedral) to the 
locations of the new particles (vertices of the tetrahedral) is equal to dl 2, where d is PVTM cell 
size). The redistribution factors are calculated using the Vorticity Redistribution Method [3.2]- 
[3.3], which is outlined below. As outlined in Section 3.1, the diffusion step is solved separately 
from the vortex convection step. Without the vortex stretching and the source terms, the 
Vorticity Transport Eqs. can be written in a non-dimensional form (Eq. 3.7), the effective 
diffusion distance, h, is defined in Eq. (3.8), where A 0 is the non-dimensional diffusion time 
step. This diffusion time step can be different from a convection time step. The redistribution 
method simulates the diffusion process within the effective diffusion distance, h, as presented in 
Eq. (3.9) by splitting a vortex particle into five smaller particles. 


d(0 1 _ 2 

— = — V 

(3.7) 

dt Re 

1 AL 


h=J- sL 

(3.8) 

V Re 

(0fx,t) = T i (j){x-x l ) 

(3.9a) 

0 ) j (x, t + AO = Y^fjTrfix - x j ) 

(3.9b) 


j=i 
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The redistribution factors, fj, are determined by using either Fourier transforms [3.2] or series 
expansion [3.3], The resulting system of equations with first order of accuracy is given in Eqs. 
(3.10), where = (x j -x i )/h, y.. = (y. - y.)/ h , and z f =(z ) ,-z,)/h . 


X /» =1 

j 

TftAj =E4iv =HfiA =° 

j j j 

=0 

j j j 

= Z-4>y =Hf'/v = 2 


(3.10) 


The system of linear equations, Eq. (3.10), is solved to obtain the redistribution factors, fy. A 
higher order of accuracy solutions can be achieved by solving additional expansion equations 
similar to Eqs. (3.10) [3.3], Once the diffusion factors are determined, each particle is split 
based on Eq. (3.11). 


a a = y Lfj° (k 


(3.11) 


3.5.2 Solving the Convection-Diffusion steps 


Since the convection and diffusion steps are solved separately, the convection and diffusion time 
steps need not necessarily to be the same. In this study, a diffusion step is performed every 20 
convection steps because the high Reynolds number makes the effective diffusion distance, h, 
very small compared to the PVTM cell size. If the diffusion step is carried out with the same 
time step as the convection step, the splitting of the particle due to diffusion effect will be 
repealed with the merging of the vortex particles (introduced to reduce the computational cost). 


3.6 Vortex stretching model for PVTM 


In addition to the reorientation of the vortex strength vector, the vortex stretching term in the 
Vorticity transport equations involves a physical stretching or thinning of the vorticity field. A 
stretching of vortex elements normally occurs when the vortex elements are subjected to a strain 
field. The direction of the stretch is governed by the invariant subspace (usually referred to as 
eigenvector) of the strain field, while the amount of stretching is governed by the principal strain 
of the strain field [3.4], In order to simulate the stretching phenomenon of a vortex element, a 
scalar variable is added to each vortex element to represent a characteristic length, L c , of the 
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element ( L c is normalized using PVTM cell size, d). Thus each vortex element has two vectors 
and one scalar quantity associated with it: (i) vortex strength vector, a, (ii) vortex location 
vector, x, and (ii) characteristic length, L c . This vortex element is subjected to a velocity field, u. 
The strain rate field, Vu, can be found by differentiating the velocity field with respect to spatial 
variables (see also Section 3.3). The principal strain of the characteristic length is found through 
Eq. 3.12: 

£ = J ma x[eig{jS7u + ±Vu T )] dt (3.12) 

Thus the equation governing the characteristic length is: 

L k c +1 =L k c { 1.0 + f) (3.13) 

where k is the time step index. The initial value of the characteristic length is d, and if the 
characteristic length of a vortex element is greater than 2 d, then the vortex element is split into 
two elements along the principle strain axis. The distance between the split elements is the 
PVTM cell size, d. After the split, the characteristic length of the split elements is reset. This 
stretching process is shown graphically in Fig. 3.3. 

3.7 Memory models for PVTM 

The implementation of the PVTM code uses either static or dynamic memory models as shown 
in Fig. 3.4. The static memory model offers fast and direct access to vortex particle information 
since the memory for all PVTM cells is allocated a priori and each cell is associated with a 
specific PVTM volume in the flow field. However, the static memory model suffers from 
inefficient usage of the memory because there are many empty cells in the PVTM domain. A 
dynamic memory model using a linked-list access is introduced to improve the efficiency of 
memory usage, but access to vorticity information in this dynamic memory model is very slow 
compared to the static memory model. The static memory PVTM model is referred to in the 
result sections as “Static PVTM”, while the dynamic memory PVTM model is referred to as 
“Dynamic PVTM”. 
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Chapter IV 


RANS/PVTM Results: Fixed-Wing Case 


In this chapter, the coupled RANS/PVTM methodology is applied to a fixed-wing case and the 
results are validated against wind tunnel test data. Comparisons with measured pressure 
distribution, loadings, and vortex parameters, and the corresponding results from the full RANS 
and coupled RANS/PVTM simulations are presented for a semi-span NACA 0015 wing. 
Although full RANS calculations are provided for reference, these full RANS calculations do not 
take advantage of higher order or refinement techniques, so they do not represent the state of the 
art in wake modeling with RANS CFD. 

4.1 Modeling Experimental Test Conditions 

A comprehensive experiment to measure the detailed progression of the tip vortex from a semi- 
span NACA 0015 wing [4.1], conducted at the NASA Ames 7- by 10-Foot Subsonic Wind 
Tunnel, is used as the reference data to validate the results from the RANS/PVTM simulations. 
The experimental setup is shown in Fig. 4.1. The experiment used a pressure-instrumented, 
untwisted semi-span NACA 0015 wing with a square tip (meaning the end is flat and 
perpendicular to the span axis). A round end cap was installed for some conditions to change 
the tip geometry of the wing from square tip to rounded tip. The experiments were conducted for 
a small range of Reynolds numbers between lxlO 6 and 3x1 0 6 . Velocity profiles across the tip 
vortex were measured at various downstream locations up to six chords behind the wing using a 
two-component laser velocimeter. The data provided by this test include (i) the chordwise 
pressure distribution along several span locations; (ii) the velocity, location, and core size of the 
tip vortex; (iii) integrated sectional lift, drag, and moment coefficients at various span locations. 

The experimental configuration and condition described in Ref. [4.1] is used to develop 
computational models of the semi-span NACA 0015 wing. The rectangular wing has a constant 
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and untwisted NACA 0015 airfoil along the span with a full-span aspect ratio of 6.6. The wing 
was mounted on a vertical supporting end plate that fitted one foot away from the side wall of the 
tunnel. This end plate prevents the formation of the trailing vortex near the mounting point and 
effectively creates an infinite wing (2D) boundary condition at the mounting point. Several 
angles of attack were tested in the experiment. For the current study, only one angle of attack of 
12° is simulated. An additional 0.51° is added to the angle of attack in the simulations as a 
correction for the closed-tunnel wall effect [4.1], The full effects of the wind tunnel walls are 
not fully accounted for in this study except near the root of the wing where only the no cross- 
flow condition is applied. The free-stream velocity corresponds to a Mach number of 0.13 and 
Reynolds number of 1 .5x1 0 6 . 

4.2 Full RANS Analysis 

The baseline or full RANS simulation results in the present study are obtained using the TURNS 
3D compressible solver [1.17], The analysis solves the RANS equations for primitive variables 
p, pV, pE inside a 3D structured grid using a finite volume approach. The grid used in the 
calculation for the square tip wing is presented in Fig. 4.2 with the dimension of 399x131x97 
grid points. The grid extends to about six chords in all directions from the wing surface, and grid 
points are distributed a priori and mostly concentrated in regions with high velocity gradient such 
as near the leading edge, boundary layer, tip vortex, and vortex sheet. The boundary condition at 
the root of the wing is set to be extrapolated from flow parameters inside the RANS domain, 
simulating the infinite wing (2D) boundary condition to match the experiment. The vorticity 
profile calculated using this full RANS analysis is shown in Fig. 4.3. It is seen that the tip vortex 
and vortex sheet are diffused very quickly before reaching the boundary of the RANS domain 
(six chords behind the wing). The rapid diffusion of the tip vortex is caused primarily by 
numerical diffusion, insufficient grid resolution (high velocity gradient regions), and the use of 
low order scheme. An improvement in the results may be obtained by using better grid 
resolution, grid adaptation scheme near high gradient region, and high order schemes. 

4.3 Coupled RANS/PVTM Analysis 

In the coupled RANS/PVTM analysis, the flow field is divided into a near-body grid and a far 
field region. The flow field in the near-body grid is resolved using the same 3D compressible 
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RANS solver and boundary conditions as described in the previous section, but the domain is 
much smaller. The RANS domain extends to about one chord length in all directions except 
behind the trailing edge, where it extends only to about a half chord length, see Fig. 4.4. The 
small extent behind the trailing edge minimizes the dissipation of vorticity to be transferred into 
the PVTM domain. This small RANS grid has a dimension of 249x131x79 grid points. The 
same 2D boundary condition at the wing root is used. In addition, the induced velocity from the 
vortex particles in the far field is included in the RANS calculation as field velocity [1.17]. The 
vorticity field resulting from this coupled RANS/PVTM calculation is presented in Fig. 4.5, and 
the vorticity isosurface of the vorticity field in the PVTM domain is shown in Fig. 4.6. These 
results demonstrate that the coupled RANS/PVTM methodology preserves the vorticity in the far 
field qualitatively quite well. 

4.4 Tip Vortex Core Parameter Identification 

The identification of the tip vortex core parameters, namely size and location, is adopted from 
the procedure presented in Ref. [4.1]. In the experiment, an approximate location of the tip 
vortex at some distance behind the trailing edge was measured using a vortex meter. Then a 
laser doppler velocimeter was used to measure the swirl velocity across the vortex core in the 
spanwise direction. The vortex core size reported is the distance between the locations of the 
maximum and minimum swirl velocity across the core. The precise location of the core was 
determined to be the mid-point between the minimum and maximum swirl velocity. 

Similar methodologies are employed for the computational results. To identify the vortex core 
parameters inside a RANS domain, the approximate location of the tip vortex core is determined 
by searching for the maximum amplitude of vorticity in a given plane behind the trailing edge. 
Once this is known, the swirl velocity across the vortex core is calculated by linear interpolation, 
and then the core size and precise location of the core are calculated in the same manner as in the 
experiment. 

A slightly different procedure is used for identifying the vortex core parameters for the PVTM 
domain, since the PVTM analysis does not calculate the velocity field directly. At the desired 
distance behind the trailing edge, properly located virtual velocity measurement planes are 
created to record the induced velocity from the vortex particle field. These measurement planes 
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are shown in Fig. 4.5 at a distance 1, 2, 4, and 6 chords behind the wing. The induced velocity 
field from the vortex particles is calculated in each measurement plane, similar to particle image 
velocimetry measurements in the experiment. The velocity fields are processed into vorticity 
with the strength vector perpendicular to the measurement plane. The vorticity field is searched 
to find the location of the maximum vorticity in order to obtain the location of the vortex core. 
The swirl velocity can be extracted directly from the measurement plane to obtain the core size 
in the same way as for the measurement and for the RANS-calculated velocity field. For the 
PVTM vortex parameter identification, 200 samples of the core size and swirl velocity field are 
averaged to obtain the reporting vortex parameters. This is necessary because the induced 
velocities from the particles are highly dependent on the proximity of the particle to the virtual 
measurement planes where the velocity is being calculated. The averaging scheme reduces the 
fluctuation of the core size and swirl velocity results due to the movement of the vortex particles 
into and out of the virtual measurement planes. 

4.5 Result comparison: Experiment, full RANS, RANS/PVTM 

Comparisons of the experimental data from Ref. [4.1] and the results obtained using the full 
RANS and the RANS/PVTM calculations are presented based on the test metric shown in Table 

4.1 . Only the static PVTM is used for the calculation of the fixed wing cases. This case metric 
is developed to provide a comprehensive representation of RANS/PVTM results that show the 
sensitivity of the RANS/PVTM results on a number of parameters including time step size (A t*), 
RANS grid boundary, PVTM cell resolution, and PVTM redistribution methods. 

Computational efficiency (time step size and number of vortex elements) is presented in Table 

4.2. All calculations are obtained using NIA high performance computing clusters. Post- 
processing of the results for plotting purposes is performed on a PC workstation using 
MATLAB. Each of the calculation presented in this report is carried out using parallel 
computing based on MPI with eight CPUs. For a comparison, the computational time and 
parameters for the full RANS calculation using four CPUs is given in Table 4.3. 

The vorticity isosurface results from the coupled RANS/PVTM calculations are presented in Fig. 
4.7 for the different cases in the test matric. The comparison of full RANS, RANS/PVTM 
results and the experimental data from Ref. [4.1] are presented in Figs 4.8-4.14, which shows 
the vorticity isosurface of the tip vortex; the chordwise pressure distribution; spanwise 
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distribution of sectional lift, drag and pitching moment; and tip vortex parameters including core 
size, vertical and spanwise locations, and swirl velocity. Overall the results from full RANS and 
RANS/PVTM calculations yield good correlation with pressure distribution data. The vortex 
parameters from the full RANS calculation do not correlate well with the data, while the vortex 
parameters from RANS/PVTM calculations correlate better with the experimental data. 

4.5.1 The effects of RANS/PVTM parameters on vorticity field 

The RANS/PVTM results in Figure 4.7 show vorticity isosurface and the location of vortex 
particles on three planes behind the wing using various combinations of RANS/PVTM 
parameters. The parameters include - (i) location of RANS grid boundary, (ii) time step size, 
(iii) PVTM cell resolution, and (iv) vortex diffusion model. The parameter variations and the 
case numbering are given in Table 4.2. It is seen from Fig. 4.7 that the effects of these 
parameters on the vorticity field in the PVTM domain qualitatively is negligible, since the results 
from all of the cases produce a very similar vorticity field behind the wing. 

4.5.2 Pressure distribution 

The comparisons of chordwise pressure distribution at four span locations are presented in Fig. 
4.8. The pressure distribution from all of the cases from RANS/PVTM calculations is similar, 
and only the results from Case 7 are shown. Overall, the pressure distribution compares well 
with the experimental data. Only the peak pressure near the leading edge is not resolved well 
from the RANS and RANS/PVTM results. The potential sources of this difference include: 
RANS grid resolution near the peak pressure, the use of low order turbulent modeling and flux 
limiter. 

4.5.3 Wing sectional loading 

Figure 4.9 shows the wing sectional loading comparison between the experimental data, full 
RANS and RANS/PVTM (Case 7) results. Only the results from Case 7 are shown since the 
differences in the RANS/PVTM results are negligible. The comparisons in sectional lift, Q, are: 
very good near the wing tip, good near the root of the wing. The sectional drag, Co, correlation 
for RANS/PVTM is very good throughout the span. The sectional moment, Cm, comparison is 
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very good near the tip, but is poor inboard, likely due to the differences in peak pressure of the 
leading edge compared to the experiment. 

4.5.4 Tip vortex core size 

The comparison of the tip vortex core size is presented in Figure 4.10. The tip vortex core size 
predicted by the RANS/PVTM method right at the release plane (at the boundary of the RANS 
grid, about half a chord behind the trailing edge) is about 30% larger than that of the test data. 
Then the core size initially increases (between 1-2 c behind the wing) and then reduces and the 
core size at six chords behind the wing is about 40% larger than the core size from the data. 
Comparing the PVTM results from all of the cases, the vortex core size from Case 5 is closest to 
the experimental data. The better correlation is probably due to the proper combination of 
PVTM parameters (small time step, smaller RANS grid in downstream direction, and the PVTM 
cell size that is similar to the data). It should be noted that this tip vortex core size is not equal to 
the PVTM cell size, for example Case 5 PVTM cell size and the resulting core size are 0.1c and 
0.14c, respectively. The core size predicted from the full RANS calculation increases rapidly 
and the core size at six chords behind the wing is more than 125% larger than the test data. This 
rapid growth in core size in the full RANS calculation is the result of numerical dissipation and 
the tip vortex moving away from the dense grid region. 

4.5.5 Tip vortex location 

The vertical and spanwise locations of the tip vortex are shown in Figs 4.11-4.12, respectively. 
Overall, the correlation of the spanwise location from RANS/PVTM is very good. The 
comparison of the vertical location is somewhat depending on the RANS/PVTM parameters. It 
is seen that the RANS/PVTM results from Case 4 give the best correlation for the tip vortex 
vertical location. It should be noted that some of the PVTM vortex vertical location results have 
sporadic movement of the tip vortex between one and two chord length behind the trailing edge. 
This is likely caused by misidentification of the vortex identification routine. For the tip vortex 
spanwise location, the results from Case 5 give the best correlation (but the vertical location 
correlation for Case 5 is not as good). 
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4.5.6 Tip vortex swirl velocity 


Comparisons of the tip vortex swirl velocity at four and six chords behind the wing are given in 
Figs. 4.13-4.14, respectively. Overall, the RANS/PVTM method predicts the swirl velocity 
outside the vortex core very well, but the comparison of the peak vortex velocity to the test data 
is fair. The peak velocity prediction is about 33% less than that of the test data at four and six 
chords behind the wing. The correlation of the velocity gradient inside vortex core at both four 
and six chords for the RANS/PVTM is better than that of the full RANS results. This is due to 
the fact that the RANS/PVTM approach preserves the vorticity better than the full RANS 
approach. 
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Chapter V 


Coupled RANS/PVTM/CSD Results: Rotary Wing Case 


In this chapter, the RANS/PVTM methodology is coupled to a CSD code to analyze a rotary 
wing case and the results are compared with hover test data. The following parameters are 
compared with the experimental data: rotor loading, vortex parameters (core size, location, 
trajectory, and swirl velocity). As a reference, a free wake hover calculation using CAMRAD II 
is included in some of the comparisons. 

5.1 Experimental test description 

An experiment to measure the detailed evolution of the rotor tip vortex in hover condition was 
conducted at the DLR Institute of Flight Systems in Braunschweig, Germany [5.1]. The test was 
performed in a limited and enclosed space (rotor preparation hall), thus the rotor was operating in 
both ground effect and recirculation flow (see Fig. 5.1). In the current study, the effects of these 
complex flow conditions are not simulated thus the rotor system is assumed to operate in an 
unconfined area. The rotor system used in this Hover Tip Vortex Structure Test (HOTIS) is the 
same as in the HART II test [5.2], Summary of the rotor and blade properties is given in Table 
5.1. The flow velocity was measured using a stereoscopic particle image velocimetry (3D-PIV) 
system. 

5.2 RANS/PVTM rotor model 

The RANS/PVTM rotor model is developed using the HART II blade properties [5.2], Four 
disconnected RANS domains are used to model the four bladed rotor system (see Fig. 2.1). Each 
RANS domain has 239x151x69 grid points (Fig. 5.2), and the grid points extend to about a half 
chord down stream and about one chord in other directions. For each RANS domain, the motion 
of each grid point is governed by the motion (three rotational and three translational 
deformations) of the quarter chord line of each blade, calculated from the CSD analysis. The 
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PVTM domain is a multi-level (six levels) domain with coarsening cell size away from the rotor 
(a representative PVTM domain with 2 levels is shown in Fig. 2.2b). The finest resolution for 
the PVTM cell is d = 0.2c, and the resolution for the PVTM cell size is defined by dx2^ LA \ 
where L is the PVTM level. 

5.3 CSD rotor model 

The structural model of the rotor is based on CAMRAD II [2.3] and is provided by NASA 
Langley researchers. The CSD model solves multi-body dynamics equations using Finite 
Element representation of the rotor system. Each blade is modeled using five structural 
elements, representing the deformation of the quarter chord line. The rotor trim is also handled 
using CAMRAD II. As a reference, a Hilly coupled trimmed calculation using CAMRAD II 
with free wake is preformed (five revolutions of the rotor wake is retained for this free wake 
calculation). 

5.4 RANS/PVTM results: untrimmed 

Using a prescribed blade motion (blade pitch angle at 0.75/? of 10°), the results from the 
RANS/PVTM (static memory model) approach is presented in Figs. 5. 3-5. 4 with three rotor 
revolutions of PVTM wake. To initialize the RANS variables, the simulation starts with only 
RANS simulation of the four-b laded rotor without releasing any PVTM wake for a quarter 
revolution (see Fig. 5.3). Without the PVTM wake and associated induced inflow, it is seen that 
the rotor thrust coefficient, Cj, increases drastically to about two times the thrust of the rotor 
with the PVTM wake. After the initial !4 revolution, the simulation starts releasing PVTM wake. 
With the induced inflow from the PVTM wake, the rotor thrust rapidly decreases over the first 
revolution and slowly increases in the second and third revolution. The rotor torque, Cq, is 
almost constant over the initial RANS startup and the first three revolutions with PVTM wake. 
The vorticity isosurface after three revolutions of the PVTM wake is shown in Fig. 5.4. It is seen 
that the startup wake (the first two revs.) is combined into a super-vortex below the rotor, and the 
strength of the tip vortex of any blade is preserved well before combining with the super-vortex. 
This shows a potential of RANS/PVTM approach to maintain the strength and compactness of 
tip vortices over one rotor revolution. 

5.5 CSD/RANS/PVTM results: coupled trimmed 
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Fully coupled trim CSD/RANS/PVTM results are presented in the section. The information 
transfer between CSD and RANS/PVTM modules occurs once every rotor revolution (loose 
coupling). As a reference, CAMRAD II (with five free wake revolutions), CII-FW, results are 
included in some of the comparisons. The following results are provided - (i) trim convergence, 
(ii) vorticity field, (iii) tip vortex swirl velocity, and (iv) tip vortex location and trajectory. 

5.5.1 Trim convergence 

Using the static PVTM model, the convergence of the trim parameter (# 75 ) is presented in Fig. 
5.5. The results include 17 resolutions of PVTM wake and involve 1 1 coupled trim loops. The 
first eight revs, of RANS/PVTM simulation is completed without retrimming to obtain a 
reasonable initial condition for the rotor wake system. The trim parameter, O 75 , from the coupled 
CSD/RANS/PVTM is about 1° higher than the DLR data and about 0.3° higher than the result 
from CII-FW. This may be because the analyses do not model the ground and recirculation 
effects. Figure 5.6 shows the variations in individual blade and rotor loading of the coupled 
results. The rotor thrust, Cj, (individual blade - Fig. 5.6a; and rotor system - Fig. 5.6b) from the 
coupled trim results converges to the thrust condition of the HOTIS test. However, the rotor 
torque, Cq, from the coupled result is 42.6% higher than the CII-FW result and is 63.2% higher 
than DLR data (CII-FW torque is 14.4% higher than DLR data). This higher torque may be 
caused by the following factors: higher 675. not modeling the ground and recirculation effects, 
and grid resolution near the blade surface to capture friction forces. 

The comparison of the sectional blade loading from the coupling iterations between CFD 
(RANS/PVTM) and CSD modules is presented in Fig. 5.7. Results from CII-FW and CII-LI 
(linear inflow model) are also included as references. It is observed that the blade sectional 
loadings (M'C/v, M 2 Cq, and M 2 Cm ) are converged after about nine loose coupling iterations. 
Comparison of the sectional normal force with CII-FW illustrates that the coupled results yield 
40-60% higher normal force near the tip of the blade (0.957?). This may be due to the differences 
in tip vortex parameter such as velocity profile and core size between the free-wake analysis and 
the PVTM approach. The blade sectional drag is much higher for the coupled results than the 
CII-FW result. Again, the factors that contribute to the higher torque value also affect the blade 
sectional drag. Blade sectional moments from the coupled results compared well with CII-FW 
result for the inboard section of the blade (r/R= 0.3-0. 7). 
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The convergence of the “delta load” during the CSD-CFD coupling iterations is given in Fig. 5.8. 
It is seen that all of the delta loadings (A M 2 , AQj M 2 , and A Cm M 2 ) are converged after about 
nine loose coupling iterations. The changes in delta loading are large in the very first coupling 
iterations and are reducing during subsequent iterations, and the changes are negligible during 
the last coupling iteration. 

5.5.2 Vorticity field 

With the coupled RANS/PVTM wake (static PVTM model) and CSD loose coupling, the 
vorticity field is shown in Fig. 5.9a. In addition to the vorticity isosurface showing the location 
of the tip vortex, Fig. 5.9a also shows virtual velocity measurement planes that rotate with the 
rotor for identifying tip vortex parameters (at 10, 20, 30,..., 150° wake ages). The tip vortex 
parameter identification procedure can be found in Appendix D. It is observed in Fig. 5.9a that 
during the first 150° wake age, the variation in the location of the tip vortex is small but the 
deviation from the mean value increases with wake age. This phenomenon is also reported in 
Ref. [5.1]. After 150° of wake age the variation in the locations of the tip vortex is noticeable 
from Fig. 5.9a. It should be noted that the tip vortex passes through RANS domains of the 
following blade and the strength of this tip vortex stays constant during and after emerging from 
the following blade RANS domains. As a comparison, the tip vortex geometry result from the 
CAMRAD II (free wake) calculation is shown in Fig. 5.9b. Overall the wake geometry from the 
coupled RANS/PVTM and CII-FW results is in good agreement in the first 270° of wake age. 
After that wake age the variation in wake geometry of the coupled RANS/PVTM wake is much 
higher than the CII-FW result. This is probably due to the differences in the solution schemes 
since free wake analysis uses non-linear iterative time stepping while the PVTM analysis uses 
explicit RK time stepping. 

5.5.3 Tip vortex swirl velocity, core size, and circulation 

Tip vortex swirl velocity is extracted from the virtual velocity measurement planes in Fig. 5.9a 
(see Appendix D for details). The tip vortex swirl velocity results (average of about 60 velocity 
measurement snapshots from the coupled RANS/PVTM wake) are presented in Fig. 5.10a for the 
wake age of 10°-150°. As the wake age increases, the swirl velocity peak and profile increases 
slightly over the range of wake age because the strength of the tip increases as it rolls in vorticity 
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from the vortex sheet. The comparison of swirl velocity from the coupled RANS/PVTM wake 
and DLR data is given in Fig. 5.10b. The data from DLR shows that the peak swirl velocity 
reduced significantly as the wake age increases from 20° -» 80° (36%) and from 80°-» 140° (38%). 
The reduction in swirl velocity is due to the effect of viscosity in and around the core of the tip 
vortex. This viscous effect is not modeled in the current PVTM methodology. At 140° wake 
age, the coupled RANS/PVTM result yields a higher swirl velocity peak than the DLR data. 
Comparison of the tip vortex circulation in Fig. 5.1 la shows that the PVTM result has larger tip 
circulation than the DLR data (possibly due to the absence of the viscosity model). The tip 
vortex core size comparison is presented in Fig. 5.11b. The predicted core size from the PVTM 
simulation is about 50% higher than the test data. It should be noted that the predicted tip vortex 
core size is different than the PVTM cell size of 0.2c. A higher resolution PVTM cell size may 
improve the prediction of the core size. 

5.5.4 Tip vortex location and trajectory 

The trajectory of the tip vortex from DLR data and coupled RANS/PVTM results are presented 
in Fig. 5.12. The coupled results include instantaneous as well as averaged locations of the tip 
vortex at 10° wake age interval. The DLR data show greater variation in the trajectory than the 
coupled RANS/PVTM results. The variation in the trajectory in the simulation is caused 
primarily by the discrete nature of the vortex particle approach that creates unsteadiness in 
velocity field, loading, and finally the location of the vortex particle. The comparison of the 
averaged tip vortex trajectory (0°-150° wake age) from the coupled results, CII-FW, and DLR 
data is shown in Fig. 5.13. It is seen that the RANS/PVTM wake trajectory is not as steep as the 
DLR data and CII-FW results, however the difference is within one chord length over the 150° 
wake age. The radial and vertical locations of the tip vortex are compared in Fig. 5.14. The 
radial location from PVTM result correlates well with that from DLR data, while the vertical 
location correlation is not as good. Overall, the differences in vertical and radial locations of the 
tip vortex between the PVTM result and DLR data is within 0.5c-0.7c. 

5.5.5 Dynamic PVTM results 

Using the converged RANS/PVTM solution from Section 5.5.1 as initial conditions, the dynamic 
PVTM model (with the vortex diffusion and stretching models) is used to simulate the rotor for 
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another two rotor revolutions. The CSD coupling (retrim) is performed every 45° to speedup the 
coupling convergence. The vorticity field from the dynamic PVTM wake is shown in Fig. 5.15. 
Comparing the vorticity field from Figs. 5.15 (dynamic PVTM with stretching model) and 5.9a 
(static PVTM without stretching model) suggests that the dynamic PVTM yields more coherent 
tip vortex almost up to one rotor revolution (a half rev. for the static PVTM). In addition the tip 
vortex is seen to be located at the center of the virtual velocity measurement planes (locations of 
these planes are set for the static PVTM results), suggesting that the change in tip vortex location 
is small. The comparison on runtime for the dynamic PVTM and static PVTM is provided in 
Table 5.2. It is seen that the dynamic PVTM run much slower than the static PVTM, due to the 
dynamic memory model that requires extra CPUs time to allocate and search the memory for 
vortex particle in the linked list. The breakdown of runtime for dynamic PVTM is presented in 
Fig. 5.16, which shows that about 60% of runtime is used for inter-CPUs communication by 
Message Passing Interface (MPI) routines (N=1000). Tip vortex trajectory results, shown in Fig. 
5.17, suggests that the dynamic PVTM with vortex stretching model improves correlation of the 
trajectory with the DLR data. This is confirmed by the correlations of the radial and vertical 
location of the tip vortex in Fig. 5.18. The radial location result track DLR data very well during 
the first 60° wake age, while the slope of the vertical location (dynamic PVTM: after 90° wake 
age) is very close to that of DLR data. However, comparison between the static and dynamic 
PVTM results (Figs 12a and Fig. 17a) shows larger scattering of the trajectory of the tip vortex 
between the wake age of 10°- 150°. This may be caused by the larger number of particles in 
dynamic PVTM producing higher level of unsteadiness in flow velocity and blade loading. For 
dynamic PVTM, the swirl velocity results are presented in Fig. 5.19. It is seen that the swirl 
velocity does decrease as the wake age increases (this trend is observed in DLR data, but not in 
the static PVTM results without the vortex stretching model in Fig. 5.10). Comparisons of the 
swirl velocity in Fig. 5.20 show good correlations of the tip vortex profile outside the vortex 
core. The core radius of the dynamic PVTM results is much larger than static PVTM results and 
DLR data. The significant reduction in peak velocity is observed in the both DLR data and 
dynamic PVTM results, but the peak velocity of the young vortex for the PVTM results (10° 
wake age, right after it is released from the RANS domain) is about half of that in DLR data (20° 
wake age). The accuracy of the PVTM result is limited by the accuracy of the RANS domain to 
preserve tip vortex strength before releasing into the PVTM domain. 
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Chapter VI 


Concluding Remarks 


This report summarizes the development and validations of a novel RANS/PVTM/CSD coupling 
framework, which is applied to simulate wake flow for fixed and rotary wing cases. The method 
employs conventional RANS solver to resolve the flow field near the bodies (rotor blades or 
wing) and gridless PVTM to provide accurate prediction of the evolution of the far field wake 
flow. For rotary wing cases, the RANS/PVTM methodology is loosely coupled to a CSD code 
that provides blade motion and trim analysis. The coupled RANS/PVTM/CSD framework 
requires proper information exchange between the modules. The results from this coupled 
framework are validated against several experimental data sets (a fixed-wing wind tunnel test 
and a hover test). 

The improvements on PVTM module were made by adding the models to simulate the effects of 
diffusion and vortex stretching. In addition, the vortex reorientation model and efficient memory 
usage model are introduced to enhance the numerical stability and efficiency of the PVTM 
module. 

For the fixed wing case, the results from RANS/PVTM method showed good correlation with 
the experimental test data. Comparisons of the pressure distribution and sectional loads were 
very good. The location of the tip vortex was also correlating well with the test data, while the 
correlations of other vortex parameters (core size and swirl velocity) were fair. 

The validation of the RANS/PVTM/CSD methodology for the rotary wing case in hover 
condition was good overall. The results showed the convergence of the coupled trim parameters 
{Cj, # 75 , C N M 2 , CcM 2 , and C m M 2 ). Rotor torque from the coupled simulation is about 50% 
higher than the test data (this was possibly caused by not modeling the ground effect and flow 
recirculating effect). The correlations of tip vortex locations were good, while that of the swirl 
velocity were fair. The tip vortex trajectory results were not as steep as the hover data but the 
difference was within one chord length over 150° wake age. The correlations on vortex 


29 



trajectory and swirl velocity are improved significantly when the vortex stretching model was 
used in the PVTM simulation. 

The accuracy of the PVTM results was seen to be highly dependent on the accuracy of the 
RANS domains to capture the location and maintain the strength of the tip vortex before 
releasing into the PVTM domain, hnprovements in the RANS domains such as higher grid 
resolution (near the tip vortex region), higher order schemes with low numerical dissipation may 
significantly improve the correlations in the future. 
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Appendix A 


Adaptive Time Step Methodology 


Adaptive time stepping for the PVTM analysis is developed to control the error associated with 
time discretization especially for particle vortex locations. The location of each vortex particle is 
obtained by solving Eq. A. 1 . 


x k+1 =x k +judt (A.l) 

where x k is the initial location and a is the induced velocity. Let us define a measure of accuracy, 

A„, , as follows (m is the numbers of time sub-step) 

A m = (X-X 0 )/(m dt), m = 1,2, M (A.2) 

An example of this measure of accuracy is shown in Fig. Al. From Fig. Al, it is seen that A m is 
an exponential function of the number of time substep, m, as follows: 

A m = A x m~ 2 (A.3) 

From Eq. A.3, the numbers of required substeps, M, to achieved a desired accuracy, £, is: 

M = yplje (A.4) 
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Appendix B 


A Two Vortex Particles Problem 


This appendix involves the analysis of a system of two vortex particles. The particles have 
strength vectors and location vectors as follows: 

aj = [ 0 1 0], xj = [ 0.05 0 0] (B.l) 

a 2 = [ 0 0 1], = [-0.05 0 0] 

Using 8 2 = (0. l) 2 /2, the induced velocity between the two particles is: 


»^2=-^ / (X2 ~' V : )X \ = [0 4.3316 0] 


(B.2) 
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The strain rate tensor, Vu, can be calculated from adding Ax, Ay, A z to the velocity evaluation 
location and re-evaluated the velocity. The strain rate tensor and its eigenvalues are: 


VUj-> 2 ~ 


0 - 43.3179 0.8661 

- 43.2942 0 0 

0 0 0 


, eig( Vu^i) = [43.3061 -43.3061 0] (B.3) 


Vu 


2->l 


0 0.8661 43.3179 

0 0 0 
43.2942 0 0 


, eig( Vu 2 ->i) = [43.3061 -43.3061 0] 


The eigenvectors Xo associated with X 0 (eigenvalue = 0) are shown in Fig. B.l. It should be noted 
that the sum of the eigenvalues is zero (continuity condition and divergence free condition). 
These eigenvectors represent the directions to which the strength of the vortex particle should be 
reoriented in order to obtain a conserved vorticity solution (see also Section 3.1.1). The direction 

of eigenvector Xo for Vu } _> ? is similar to CCj, and the eigenvector Xo for Vu 2 -> i is pointing in the 
direction of a 2 . This is true for any pair of vortex particles. Thus in general, the reorientation of 
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the vortex strength vector should be in the direction of the eigenvector, Xo- To simplify the 
process for vorticity field with more than one particles, the vortex stretching term is modified 
from Eq. B.4 (direct differentiation of Eq. B.2) to B.5 to take into account the proper 
reorientation of each vortex particle. Equation B.6 is introduced to adjust the reorientation speed 
of the strength vectors. 


a p . Vu _y 1 a r *a t | 3(a f • (x f - x, ) x or, )(x f - x,. ) 

j47T(\ P 1 2 f\ n 12 

(J X 


a m Vu i ='Y J 
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Using Eqs. (B.4-B.6), the vortex trajectory of this two vortex particle system is shown in Fig. 
B.2. It should be noted that the two vortex particles realign themselves properly after several 
time steps using Eqs (B.5 and B.6), and the total vorticity is conserved (\ai N \ = \(Xi°\, \o( 2 N \ = |C6°|). 
The distance between the two vortex particles is presented in Fig. B.3. It is seen that the distance 
between the vortex particles is increasing unboundedly when Eq. B.5 and B.6 are used. To 
mitigate this unbounded drifting (associated with explicit time integration scheme), an attraction 
term is introduced in the velocity calculation (Eq. B.2) as follows: 


u(Y)= Z" 


4/r 


(Y -x,)xa t + 0.05 (Y - x, )|(F-x,)x a, 
\y-x, \ 2 +S 2 f 


(B.7) 


Using Eqs. (B.5) and (B.7), the distance between the two vortex particles is similar to that 
obtained using Eq. (B.4), see Fig. B.3. 


It is recognized that Eq. (B.7) is not derived from any physical phenomena, and the following 
section serves to identify the source(s) of the modification introduced in Eq. (B.7). One possible 
source is the accumulated error from the time integration routine (explicit Rouge-Kutta fourth 
order). As a comparison the following time integration routines are used and the results are 
compared in Fig. B.4 (using Eq. (B.2) and (B.6), and a different initial vortex strength and 
location from the previous section). 

(a) Explicit Rouge-Kutta (6 th order) 

(b) Implicit (preditor-corrector) Adam-Bashford (2 ud , 3 rd , and 4 th order) 
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The comparison shows no improvement using the higher order explicit integration schemes, but 
the introduction of the implicit time integration routines can actually improve the “particle 
drifting” problem. With these results, it can be concluded that the “particle drifting” problem 
can be solved using higher order implicit integration schemes with a sufficiently small time step. 
However, in practice using the implicit time integration routine increases the memory 
requirement by a factor of two, three, or four depending on the required order of accuracy. This 
requirement is impractical when the system involves very large number of particle (100000+). 
In addition, the time step required is extremely small making the simulation very 
computationally expensive. 

Thus, for engineering purposes, a modification is introduced to Eq. (B.2) to include the attraction 
term as shown in Eq. (B.8) as a way to reduce the discretization error from the explicit time 
integration routine used (RK4). This attraction term is used only for calculating velocity from 
particle to another particle, and not in the induced velocity calculation for coupling with the 
RANS domains. 


»(T) = £ --T- 


(Y 


*.) XQr . +7fE%f|( y -*/)xg,- 


(jy-x ,.| 2 + s 2 f 


(B.8) 


where the attraction coefficient, rj, is a variable controlling the time discretization error. It is 
assumed that rj is a function of the time step. At, and the optimal values of r) can be determined 
by running the simulation adjusting IJ such that the “particle drifting” is negligible (with RK4 
scheme for various At). As a function of At, the optimal attraction coefficients, tji and t) 2 , are 
shown in Fig. B.5. A trajectory of a two particles system with and without the attraction term is 
shown in Fig. B.6. It is observed that the inclusion of the attraction term, rji, significantly 
reduced the “particle drifting” problem, however the distance between the particles is still 
increasing. A small offset is added to tji resulting in 1)2 to make the “particle drifting” problem 
negligible (see Figs. B.5- B.7). 
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Appendix C 


Single Bladed Rotor Cases 


This appendix presents results for single bladed rotor cases that demonstrate the effects of some 
PVTM modeling methods and parameters including (i) dynamic PVTM methodology, (ii) PVTM 
stretching model, and (iii) high resolution PVTM. 

C.l Dynamic PVTM 

The single-bladed hover results using the static PVTM and dynamic PVTM methodology are 
shown in Figs. C.1-C.2. The vorticity fields presented in Figs. C.l show that the static PVTM 
and dynamic PVTM results are equivalent. The rotor loading comparison in Fig. C.2 confirms 
that the rotor loading is the same using the two approaches. The effect of the attraction term 
(Appendix B) is presented in Fig. C.3. The tip vortex from PVTM wake with the attraction term 
is much more compact and coherent than the PVTM result without the attraction term. 

C.2 PVTM stretching model 

With the vortex stretching methodology (Section 3.6), the vortex particle releases from a single- 
bladed rotor is shown in Fig. C.4 (color representing the characteristic length, L c ). The figure 
shows the characteristic length of the vortex sheet at 2 azimuthal positions. It is observed that 
the vortex particles near the tip vortex undergo splitting due to vortex stretching very quickly 
after the release (about 5°-10° wake age) the rest of the vortex sheet reaches the threshold for 
splitting after approximately 20°-45° wake age. The splitting of the vortex particle in the vortex 
sheet continues after the initial splitting. Comparison of the particle field with and without 
vortex stretching model is presented in Fig. C.5. The resulting particle field with the vortex 
stretching model is seen to be much more uniform, especially in the vortex sheet near the tip 
vortex. Comparisons of rotor loading seen in Fig. C.6 shows that the PVTM results with the 
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vortex stretching model yield a similar rotor loading without the stretching model. 
Computational efficiency of the PVTM with vortex stretching model is summarized in Table 
C. 1 . With the vortex stretching model, both the numbers of particles and the runtime increase by 
about 300% from the baseline on stretching case. 

C.3 High resolution PVTM wake 

With the development of the dynamic PVTM methodology using the dynamic memory model, a 
simulation with high resolution PVTM wake is possible. The differences in particles released in 
one time step at various PVTM resolutions can be seen in Figs. C7-C.8. By observing the 
clustering of the vortex particles near the tip vortex, it is seen the high resolution PVTM wake 
result in Fig. C.8b provides better tip vortex profiles (location, strength distribution) than the low 
resolution PVTM wake. The high resolution PVTM wake also affects the resolution of the 
coupling methodology between the RANS and PVTM domains. For example, coupling a RANS 
domain to PVTM domain (or other RANS domains) involves a conversion from the RANS 
vorticity field (Fig. C.9) into a particle representation of the field (Fig. C.10-C.11). A high 
resolution PVTM conversion (Fig. C.ll) yields four times as many particles than a lower 
resolution PVTM conversion (Fig. C. 10), but again provides better vorticity profile. Comparison 
of vorticity field, rotor loading, and computational efficiency of various PVTM wake resolution 
after one rotor revolution are presented in Figs. C.12 and C.13, and Table C.2, respectively. 
Overall the high resolution PVTM wake improves the vorticity profile (without changing the 
rotor loading), while increasing the numbers of vortex particles and the runtime. 
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Appendix D 


Rotor Tip Vortex Parameter Identification 


This appendix outlines the procedure to identify the tip vortex parameters, i.e. core size, swirl 
velocity, and location. Starting with the 3D induced velocity recorded in a measurement plane, 
the 3D velocity is projected onto the measurement plane to yield a 2D velocity field, and the 2D 
velocity is processed to calculate the vorticity perpendicular to the measurement plane as shown 
in Fig. D.l. The vorticity is calculated using a fourth order central difference scheme, and a 2-D 
gaussian filter is applied to the processed vorticity results to reduce the numerical noise. An 
estimate of the location of the tip vortex, X t , p 2 , can be found by searching for the maximum 
vorticity in the measurement plane. Using this estimated location as a starting point, a more 
sophisticated scheme is used to obtain more accurate location of the tip vortex, X t i p ', using the 
velocity profile search. The scheme identifies the “I index” of the tip vortex by searching the 
vertical velocity, Vi, for the mid point between minimum and maximum velocity and its 
derivative, dVi/dl, for the maximum value of the derivative as shown in Fig D.2. Similarly, the 
“J index” of the tip vortex is identified by searching the horizontal velocity, Uj, and its 
derivative, dUj/dJ. The final location of the tip vortex, X t i P , is simply the average between Xtip 1 
and X t i p 2 . Once the location of the tip vortex is identified, the horizontal (U) and vertical (V) 
swirl velocities are obtained using ID interpolation in the appropriate direction. The final swirl 
velocity, V 9 , is the average between the horizontal and vertical swirl velocities (Fig. D.3). 
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Tables 


Table 4.1 : RANS/PVTM case metric 


Cases 

A t* 

CFD grid boundary 

PVTM cell 
resolution 

PVTM redistribution 
methods 

Wrap 

around 

direction 

(chord) 

Down 

stream 

direction 

(chord) 

Level 0 
(chord) 

Level 1 
(chord) 

Stretching 

model 

Diffusion 

model 

0 

0.02 

1 

0.5 

0.1 

0.4 



1 

0.02 

1 

0.5 

0.05 

0.2 



2 

0.02 

1 

0.5 

0.05 

0.2 


X 

3 

0.02 

1 

1 

0.05 

0.2 



4 

0.02 

1 

1 

0.05 

0.2 


X 

5 

0.01 

1 

0.5 

0.1 

bei 





1 

0.5 

0.05 




7 

0.01 

(adaptive) 

1 

0.5 

0.05 

0.2 




non-dimensinalized 


by c/a, (c: wing chord, a: speed of sound) 


Table 4.2: Number of vortex particles and calculation time for the RANS/PVTM 
calculations using 8 CPUs 


Cases 

CFD-PVTM 

Total 

Time Steps 

PVTM vortex 
particles 

Calculation time of the last time step 

Level 0 

Level 1 

CFD (%) 

PVTM (%) 

Total (sec.) 

0 

9000 

1355 

76 

35.02 

64.98 

98.81 

1 

9000 

5415 

291 

46.50 

53.50 

111.96 

2 

9000 

5484 

286 

46.48 

53.52 

112.11 

3 

9000 

4159 

307 

47.27 

52.73 

■BOH 

4 

9000 

4171 

295 

47.02 

52.98 


5 

17000 

1712 

71 

34.34 

65.66 

119.56’ 

6 

17000 

4978 

290 

31.51 

68.49 

87.54 

7 

17000 

5046 

296 

33.51 

66.49 

85.19 


*A dillerent computing cluster is used, and tins cluster is running about 20% slower 


Table 4.3: Full CFD case parameters and calculation time using 4 CPUs 





CFD grid boundary 


Cases 

CFD Total 
Time Steps 

A t* 

Wrap around 
direction 
(chord) 

Down stream 
direction 
(chord) 

Calculation time of the 
last time step 

0A 

10000 

0.02 

6 

6 

73.78’ 


*A different computing cluster is used. and tins cluster is miming about 20% slower 
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Table 5.1: Properties of the HOTIS (HART II) model rotor 


Property 

Value 

Number of blades 

4 

Radius 

2.0m 

Root cutout 

0.44m 

Chord 

0.121m 

Solidity 

0.077 

Airfoil 

NACA230 1 2mod 

Linear twist 

-8.0deg/R 

Nominal speed 

109 rad/s 

Tip mach number 

0.633 


Table 5.2: PVTM runtime com] 

parison 

Model 

Number of 
particles 

Number 
of CPUs 

Runtime 
RANS (%) 

Runtime 
PVTM (%) 

Runtime per time 
step (sec) 

Static 

PVTM 

101649 

12 

42.15 

57.85 

105.75 

Dynamic 

PVTM 

84578 

8 

4.61 

95.39 

1301.72 


Table C.l: Comparisons of computational efficiency of static and dynamic PVTM (with 
and without vortex stretching) 


Case 

d/c 

Particles 

Run time (sec) 
Last time step 

CPUs 

Static PVTM 

0.20 

8183 

66.27 

4 

Dynamic PVTM 
(no stretching) 

0.15 

16481 

92.93 

8 

Dynamic PVTM 
(with stretching) 

0.15 

45585 

330.55 

8 


Table C.2: Comparisons of computational efficiency of static and high resolution PVTM 


Case 

d/c 

Particles 

Run time (sec) 
Last time step 

CPUs 

Static PVTM 

0.20 

8183 

66.27 

4 

Dynamic PVTM 

0.20 

11411 

59.19 

4 

0.15 

16481 

92.93 

8 

0.10 

29018 

336.68 

8 
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Figures 


Computational Fluid Dynamics 


Computational Structural Dynamics 

(CFD) 


(CSD) 



Figure 1.1: Overview of coupled PVTM/RANS/CSD methodology 
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(a) fixed-wing case 




(b) rotor case 


Figure 2.1: RANS grid for near filed flow calculation: 
(a) fixed-wing case and (b) rotor case 
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Figure 2.2: Two dimensional view of representative 3D multi-level PVTM domains: 

(a) fixed-wing case and (b) rotor case 



Figure 2.3: (a) Particle representation of a continuous vorticity field, (b) Merging of vortex 
particles in the same PVTM computational cell 
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(a) released particles (convection) 



(b) released particles (grid movement) 


VTM cell 

(c) released particle after merging 

Figure 2.4: Vortex particle released from RANS domain: (a) from convection process, 
(b) from grid movement, and (c) after merging 



Figure 2.5: Induced velocity from PVTM domain on interpolation grid before using 3D 
linear interpolation to calculate induced inflow on RANS grid 
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Figure 2.6: Loose coupled trim procedure for coupling RANS/PVTM to CAMRAD II 


Vortex particle 



Figure 3.1 : Interpolation basis points for calculating induced velocity within RK scheme 

and local strain tensor 
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Figure 3.2: Vortex redistribution model than simulate diffusion process 




Principal strain field 

(a) t = t 0 

Figure 3.3: Vortex stretching from invariant strain tensor after k time steps 



(a) static memory model (b) dynamic memory model 

(Static PVTM) (Dynamic PVTM) 

Figure 3.4: Memory map for PVMT implementation with (a) static and (b) dynamic 

memory models 
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Figure 4.1 : Experimental set up for the NACA 0015 wing in wind tunnel tests 

(from Ref. [4.1]) 



Figure 4.2:Grid for full RANS calculation of the square tip wing (5M grid points) 
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Figure 4.3: Vorticity profile calculated using full RANS calculation for square tip wing 



Figure 4.4: RANS grid for coupled RANS/PVTM calculation of the square tip wing (2.5M 

grid points). 
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(b) 

Figure 4.5: Vorticity profile calculated using coupled RANS/PVTM method showing 
velocity planes and vortex particles in the PVTM domain. 
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Figure 4.6: Vorticity isosurface and particle trace on three vertical planes of the wake 
behind a square tip wing simulated using the coupled RANS/PVTM method 
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Figure 4.7: Vorticity isosurface and particle trace on three vertical planes of the wake 
behind a wing simulated using RANS/PVTM method 
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Pressure Coefficient Pressure Coefficient 




(a) 0.974R 


(b) 0.944R 




(c) 0.845R 


(d) 0.773R 


Figure 4.8: Airfoil surface pressure distributions near the wing tip (case 7) 
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Figure 4.9: Sectional lift, drag, and moment coefficients along wing span (case 7) 
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Tip vortex core size (c) Tip vortex core size (c) Tip vortex core size (c) Tip vortex core size (c) 










Distance behind wing trailing edge (c) Distance behind wing trailing edge (c) 


Figure 4. 10: Comparison of tip vortex core size 
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Distance behind wing trailing edge (c) 


Distance behind wing trailing edge (c) 


Figure 4. 1 1 : Comparison of tip vortex core vertical location 
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Distance behind wing trailing edge (c) Distance behind wing trailing edge (c) 


Figure 4.12: Comparison of tip vortex core spanwise location 
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Distance from vortex core (mm) Distance from vortex core (mm) 


Figure 4. 13: Comparison of tip vortex swirl velocity (V/V^) at 4 chords behind the wing 
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Distance from vortex core (mm) 


Distance from vortex core (mm) 


Figure 4. 14: Comparison of tip vortex swirl velocity (V/V^) at 6 chords behind the wing 


59 



Figure 5.1: HOTIS test configuration (from Ref. [5.1]) 
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Figure 5.2 RANS grid for one of HOTIS blades 



Figure 5.3: Rotor loading (untrimmed HOTIS, 3 revs) 
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Figure 5.4: Vorticity isosurface of the PVTM vorticity field (untrimmed HOTIS, 3 revs) 
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Figure 5.5: Variation of trim parameter, 675, for loose coupling results for HOTIS rotor 
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(b) Rotor system 

Figure 5.6: Variation in CFD loadings (Cr and Cq) for HOTIS rotor in hover 
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x 10' 3 







Figure 5.7: Comparisons of blade loading from CFD and CSD analyses during 9 trim 
iterations [Free wake (FW) and linear inflow (LI) values are added for completeness] 
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0 3 04 0 5 0.5 07 0 8 0 9 1 



Radial station (r/R) Radial station (r/R) 



Figure 5.8: Comparisons of “Delta loading” during nine trim iterations 

(e.g. ACn = Cn CFD - Cn LI ) 
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Virtual velocity measurement plane 


Projection of blade tip 
on a measurement plane 


4 

t 


(a) Coupled RANS/PVTM/CSD result 



(b) CAMRAD II (free wake) result 


Figure 5.9: Vorticity field from the converged HOTIS results 
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Swirl velocity (Vq/QR) 



(a) PVTM results 



(b) Comparison with DLR data 


Figure 5. 10: Tip vortex swirl velocity: (a) converged RANS/PVTM/CSD results and (b) 

comparison with DLR data 
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Circulation (I7QRc) 



0r075 


0.15 


Ot125 


O Circulation (psi = 20 deg) DLR 
— ♦ — Circulation (psi = 80 deg) DLR 
■ Circulation (psi = 140 deg) DLR 
Circulation (psi = 140 deg) PVTM 


-0.3 -0.2 -0.1 0 0.1 0.2 0.3 

Distance from vortex core (c) 

(a) tip vortex circulation (r = JVedl = JVer d0) 



Figure 5.11: Comparisons of (a) tip vortex circulation and (b) tip vortex core size 
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2 0.4' 1 4 j 


a 



Radial distance from tip (c) 
(b) DLR data (from Ref. [5.1]) 


Figure 5.12: Tip vortex trajectory from (a) RANS/PVTM/CSD result and (b) DLR data 
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Vertical location (c) 


0.5 



—■—Vortex trajectory (DLR) 
—♦—Vortex trajectory (PVTM) 
—^•Vortex trajectory (CII-FW) 


Radial location from blade tip (c) 


Figure 5.13: Comparison of tip vortex trajectory (0°-150° wake age, O: indicates the first 

blade passage of 90°) 
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Vertical location from tip (c) Radial location from tip (c) 
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Wake Age (deg.) 

(a) Radial location from tip (c) 



0.5 


0 


- 0.5 


-1 



Wake Age (deg.) 

(b) Vertical location from tip (c) 


Figure 5.14: Comparison of tip vortex location (0°-150° wake age) 
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Runtime per particle (sec/N) 



(a) Dynamic PVTM (b) CAMRAD II (free wake) 

Figure 5.15: Comparison of vorticity field: dynamic PVTM and CII-FW 
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Figure 5.16: Runtime for dynamic PVTM 
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Vertical location (c) 



(a) Dynamic PVTM 



Vortex trajectory (DLR) 

♦ Vortex trajectory (Static PVTM) 
—A — Vortex trajectory (Dynamic PVTM) 
— * — Vortex trajectory (CII-FW) 


2.5 


- 0.5 


- 1.5 


Radial location from blade tip (c) 


(b) Comparison with other results and data (O: indicates the first blade passage of 90°) 
Figure 5.17: Tip vortex trajectory from RANS/PVTM/CSD result (dynamic PVTM) 


73 


Vertical location from tip (c) Radial location from tip (c) 


0.5 



Wake Age (deg.) 

(a) Radial location from tip (c) 



Figure 5.18: Comparison of tip vortex location (Dynamic PVTM) 
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Swirl velocity (Ve/QR) 


0.15 



Figure 5. 19: Tip vortex swirl velocity from converged RANS/PVTM/CSD result 

(Dynamic PVTM) 



Figure 5.20: Comparisons of the tip vortex swirl velocity (Dynamic PVTM) 
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2.00E-05 



Number of integration substep, m 

Figure A1 : Measure of accuracy, A m , as a function of number of sub-steps, m 




(a) Vu 2 ->i (b) Vuj +2 

Figure B.l : Eigenvectors, %o , associated with Ao for the two particles system 
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Distance between particles 




Figure B.2: Trajectory for a system of two vortex particles using the stretching term from: 

(a) Eq. B.4, (b) Eq. B.5, and (c) Eq. B.6 



Figure B.3: Distance between the two vortex particles with various stretching model 
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Distance between particles 



Figure B.4: Comparison of results using explicit and implicit time integration routines 



Time step (At) 

Figure B.5: Optimal attraction coefficient as a function of the time step 
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Distance between particles 





(a) RK6 (b) T|i (c)Tl 2 

Figure B.6: Trajectory of two vortex particle system: (a) RK6-no attraction term. 


(b) RK4-r|i ,(c) RK4-T12 (A/ = 0.5 tt/180) 
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(a) Static PVTM (8178 particles) 



(b) Dynamic PVTM (11161 particles) 

Figure C.l: Comparison of vorticity field after 1 rev.: (a) static PVTM and (b) dynamic 

PVTM (d = 0.2c) 
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(a) with attraction term (b) without attraction term 

Figure C.3: Comparisons of vorticity field from PVTM with and without attraction term 
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(a) with vortex stretching (b) without vortex stretching 

Figure C.5: Evolution of vortex sheet and tip vortex shed from a hovering rotor with and 

without vortex stretching ( d = 0.15c) 
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16 - 



Figure C.6: Comparisons of the blade loading with and without stretching model 



(a) Static PVTM (103 particles) (b) Dynamic PVTM (93 particles) 

Figure C.7: Particles released at one time step from HOTIS rotor in hover ( d = 0.2c) 





(a) d = 0.1c (191 particles) (b) d = 0.05c (306 particles) 

Figure C.8: Particles released at one time step from HOTIS rotor in hover 

(dynamic PVTM) 
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Figure C.9: Vorticity field in RANS domain before conversion 




(a) Static PVTM (1001 particles) (b) Dynamic PVTM (902 particles) 

Figure C. 10: Particle vortex representation of the vorticity field after conversion ( d = 0.2c) 



Dynamic PVTM (3999 particles) 

Figure C. 1 1 : Particle vortex representation of the vorticity field after conversion (d = 0. 1 c) 
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(a) Static PVTM (8178 particles, d = 0.2c) 



(c) Dynamic PVTM (16481 particles, d = 0.15c) 



(d) Dynamic PVTM (29018 particles, d = 0.10c) 

Figure C.12: Vortex particle field after 1 rotor revolution using static and dynamic PVTM 

implementations 
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Figure D. 1 : Vorticity on the measurement plane with the wake age of 10° 
(Vorticity strength is perpendicular to the measurement plane) 
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Swirl velocity 






(a) (b) 

Figure D.2: Identification of tip vortex location using (a) vertical swirl velocity and (b) 

horizontal swirl velocity 



Figure D.3: Swirl velocity of the identified tip vortex 
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